home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / Libraries / DCLAP 6d / dclap6d / SeqPups / appsrc / autoseq.src / CTraceFile.H < prev    next >
Text File  |  1996-07-05  |  10KB  |  397 lines
  1. //    ============================================================================
  2. //    CTraceFile.H                                                    80 columns
  3. //    Reece Hart    (reece@ibc.wustl.edu)                                tab=4 spaces
  4. //    Washington University School of Medicine, St. Louis, Missouri
  5. //    This source is hereby released to the public domain.  Bug reports, code
  6. //    contributions, and suggestions are appreciated (to email address above).
  7. //    -    -    -    -    -    -    -    -    -    -    -    -    -    -    -    -
  8. //    CTraceFile represents chromatogram data from automated DNA sequencers.  It
  9. //    currently supports reading and writing files in ABI and SCF formats.
  10. //    Thanks to LaDeana Hillier for providing the code from which these file
  11. //    formats were inferred.
  12. //
  13. //    NOTES
  14. //    !    This source makes no attempt to accommodate machine independent I/O
  15. //        because it was more important to tackle other issues first.  Nonetheless
  16. //        I have tried to implement most I/O functions in a way that will make
  17. //        migration to machine-independent I/O.
  18. //    ?    bottom flag is unimplemented.  the seqIO implementation of reading 
  19. //        each trace in reverse for bottoms is, ahem, interesting.  I'd suggest 
  20. //        using the template function Invert in RInlines to invert both sequences
  21. //        and traces.
  22. //    *    CTraceFile is generic in that it can represent several file formats.
  23. //        Some fields of the abi format have been omitted and are copied into the
  24. //        comments field during reading.  Before this class can correctly write
  25. //        abi files, the field integrity must be maintained.  This may be done
  26. //        by (1) making a derived class which declares these members and
  27. //        appropriate read & write methods, or (2) make this class a superset of
  28. //        all file formats.
  29. //    *    SCF doesn't currently read or write comments
  30. //
  31. //    MODIFICATION HISTORY
  32. //    93.11.11    Reece    First release
  33. //    ========================================|===================================
  34.  
  35. #ifndef _H_CTraceFile                        // include this file only once
  36. #define _H_CTraceFile
  37.  
  38. #include    <stddef.h>
  39. //#include    "CTrace.H"
  40. #include    "CTrace.C"  // dgg
  41. #include    "CPeakList.H"
  42. #include    "DNA.H"
  43. #include    "FileFormat.H"
  44. #include    "RInclude.H"
  45.  
  46. typedef        double xform_mtx[NUM_BASES][NUM_BASES]; // See Transform(), below
  47.  
  48. static const double MAX_SCF_TRACE_VALUE = 255.0;
  49.  
  50. class    CTraceFile
  51.     {
  52.     public:
  53.     enum    strand_t
  54.         { top, bottom };
  55.  
  56.     enum    error_t
  57.         {
  58.         noError,                            // no error occurred
  59.         memError,                            // memory couldn't be allocated
  60.         ioError,                            // file/strm access error
  61.         fileExistsError    ,                    // file already exists
  62.         fileDoesntExistError,                // no such file
  63.         emptyError,                            // I'm empty and can't do that
  64.         fmtError,                            // I don't understand data format
  65.         unkFmtError,                        // can't figure out the format
  66.         fmtNotSuppError,                    // format not yet supported
  67.         peakPickError                        // error in peakpicking
  68.         };
  69.  
  70.     //
  71.     // Instance variables
  72.     //
  73.     private:
  74.     static vrsn    version;
  75.     error_t    error;                            // error flag
  76.  
  77.     char*    filename;                        // filename
  78.     format_t    nativeFormat;                // format used to read this data
  79.     strand_t    whichStrand;                // which strand are we looking at?
  80.  
  81.     CTrace<trace_t>* trace[4];                // the traces
  82.     size_t        numPoints;                    // # of points in traces
  83.     tracepos    leftCutoff;                    // left cutoff
  84.     tracepos    rightCutoff;                // right cutoff
  85.     tracepos    primerPosition;                // start of primer
  86.     CPeakList    peaks;                        // assimilated (ACGT) list of peaks
  87.     trace_t        min;                        // min,
  88.     trace_t        max;                        // max,
  89.     trace_t        mean;                        // & mean for the 4 traces
  90.  
  91.     size_t        numBases;                    // # of bases
  92.     char*        sequence;                    // the bases as called by mfr
  93.     tracepos*    basePositions;                // base # <-> trace point # rel'n
  94.  
  95.     char*        comments;                    // miscellaneous comments
  96.  
  97.  
  98.     //
  99.     // Methods
  100.     //
  101.     public:
  102.     vrsn&        Version();                    // return class version
  103.     error_t        Error(void);
  104.     void        Error(error_t new_error);    //   set/clear error
  105.  
  106.                 CTraceFile(                    // constructor
  107.                     const char* fn=NULL,
  108.                     format_t fmt=unknown);
  109.                 ~CTraceFile(void);            // destructor
  110.  
  111.     error_t        Allocate(                    // allocate traces,
  112.                     size_t points,            // base positions,
  113.                     size_t bases,            // and comment string
  114.                     size_t comment);
  115.     void        Deallocate(void);            // deallocate space
  116.  
  117.     void        NumPoints(size_t np);        // set and...
  118.     size_t        NumPoints(void);            //   get # of points in traces
  119.     trace_t        Mean(void);                    // get mean,
  120.     trace_t        Min(void);                    //   min, &
  121.     trace_t        Max(void);                    //   max values over all 4 traces
  122.     void        CalculateStats(void);        // calc. min/max/mean for all traces
  123.  
  124.     CPeakList&    Peaks(void);                // get the list of peaks
  125.     size_t        NumPeaks(void);                // # of peaks picked 
  126.     void        NumBases(size_t nb);        // set and...
  127.     size_t        NumBases(void);                //   get the number of bases
  128.     void        Sequence(char* seq);        // set and...
  129.     char*        Sequence(void);                //   get the sequence
  130.     void        Comments(char* newComment);    // set and...
  131.     char*        Comments(void);                //   get the comment
  132.  
  133.     CTrace<trace_t>*
  134.                 operator[](enum_t whichTrace);    // trace selector
  135.  
  136.     tracepos    LeftCutoff(void);            // get and...
  137.     void        LeftCutoff(tracepos t);        //   set left cutoff
  138.     tracepos    RightCutoff(void);            // get and...
  139.     void        RightCutoff(tracepos t);    //   set right cutoff
  140.  
  141.     tracepos    PrimerPosition(void);        // return the position of the
  142.                                             // primer (if stored in file)
  143.     tracepos*    BasePositions(void);        // return base positions array
  144.  
  145.     error_t        Transform(xform_mtx& matrix);
  146.                 // Applies a 4x4 transformation/orthogonalization matrix to the
  147.                 // 4 traces, producing a new set of traces which replaces the
  148.                 // existing set.
  149.                 // Transformation matrics are expected to be 4x4 matrix of the
  150.                 // form:
  151.                 //     M =         [,A]    [,C]    [,G]    [,T]
  152.                 //         [A,]    mAA     mAC     mAG     mAT
  153.                 //         [C,]    mCA     mCC     mCG     mCT
  154.                 //         [G,]    mGA     mGC     mGG     mGT
  155.                 //         [T,]    mTA     mTC     mTG     mTT
  156.                 // The general equation for the resulting traces is:
  157.                 // R = M O  <==>  R(T,i) =      sum     [ mTS x O(S,i) ]
  158.                 //                          S in {ACGT}
  159.                 // where R is the resulting vector of 4 traces
  160.                 //       O is the original vector of 4 traces
  161.                 //       M is the 4x4 ({ACGT}x{ACGT}) matrix whose elements
  162.                 //         m(i,j) are the cross-term contributions of channel j
  163.                 //         to channel i
  164.                 //       S & T are trace identifiers (Source & Target)
  165.                 //         in {A,C,G,T}
  166.                 //       i loops over the indices of the trace
  167.  
  168.  
  169.     error_t        PickPeaks(
  170.                     double PeakMeanCoefficient = 1.5,
  171.                     double ZeroThreshold = 0.0);
  172.                 // Calls PickPeaks on each trace in the tracefile, then
  173.                 // assimilates the peaks into a list of peaks for the entire
  174.                 // set sorted by index in trace.  No pruning in performed.
  175.  
  176.     error_t        AssimilatePeaks(void);
  177.                 // Performs a merge sort of the list of peaks stored in each
  178.                 // trace's peak list and puts the result in this tracefile's
  179.                 // own (assimilated) peak list.
  180.  
  181.     error_t        PrunePeaks(
  182.                     tracepos minSeparation,
  183.                     ostream* os=NULL);
  184.                 // Prunes peaks by doing a pairwise comparison of adjacent peaks
  185.                 // and discarding the less-probably peak of each pair which is
  186.                 // separated by less than minSeparation.  If os is specified,
  187.                 // a list of pairwise comparisons which resulted in the
  188.                 // removal of a peak is output.
  189.  
  190.     //
  191.     // I/O methods
  192.     //
  193.     error_t        Read(
  194.                     const char* filename,
  195.                     format_t fmt=unknown);
  196.                 // Top-level read routine.  Attempts to read the named file
  197.                 // in the specified format, or in the format guessed by
  198.                 // FileFormat() (see FileFormat.H).
  199.     error_t        Write(
  200.                     const char* filename,
  201.                     format_t fmt=unknown);
  202.                 // Top-level write routine.  Attempts to write the named file
  203.                 // in the specified format, or in the native format if
  204.                 // unspecified.
  205.  
  206.  
  207.     private:
  208.                 // The following are private I/O methods.  Users should not
  209.                 // need access to any of these.
  210.     error_t        Read(                        // read from file, format known
  211.                     FILE* fp,
  212.                     format_t fmt);
  213.     error_t        Write(                        // write to file, format known
  214.                     FILE* fp,
  215.                     format_t fmt);
  216.     error_t        ReadSCF(FILE* fp);            // Read as SCF file
  217.     error_t        WriteSCF(FILE* fp);            // Write as SCF file
  218.     error_t        ReadABI(FILE* fp, short whichSet=0);    // Read as ABI file
  219.     error_t        WriteABI(FILE* fp);            // Write as ABI file
  220.     error_t        ReadALF(FILE* fp);            // Read as ALF file
  221.     error_t        WriteALF(FILE* fp);            // Write as ALF file
  222.  
  223.     public:
  224.     friend ostream& operator<<(ostream& os, CTraceFile& ctf);
  225.                 // Places an easy-to-read summary of the tracefile contents
  226.                 // on the specified output stream.
  227.     };
  228.  
  229.  
  230. //
  231. //    INLINE FUNCTION DEFINITIONS
  232. //    ===========================
  233. //
  234.  
  235. inline
  236. vrsn&
  237. CTraceFile::Version()
  238.     {
  239.     return version;
  240.     }
  241.  
  242. inline
  243. void
  244. CTraceFile::Error(error_t new_error)
  245.     {
  246.     error = new_error;
  247.     }
  248.  
  249. inline
  250. CTraceFile::error_t
  251. CTraceFile::Error(void)
  252.     {
  253.     return error;
  254.     }
  255.  
  256. inline
  257. void
  258. CTraceFile::NumPoints(size_t np)
  259.     {
  260.     numPoints = np;
  261.     }
  262. inline
  263. size_t
  264. CTraceFile::NumPoints(void)
  265.     {
  266.     return numPoints;
  267.     }
  268. inline
  269. trace_t
  270. CTraceFile::Mean(void)
  271.     {
  272.     return mean;
  273.     }
  274. inline
  275. trace_t
  276. CTraceFile::Min(void)
  277.     {
  278.     return max;
  279.     }
  280.  
  281. inline
  282. trace_t
  283. CTraceFile::Max(void)
  284.     {
  285.     return max;
  286.     }
  287.  
  288. inline
  289. tracepos
  290. CTraceFile::LeftCutoff(void)
  291.     {
  292.     return leftCutoff;
  293.     }
  294. inline
  295. void
  296. CTraceFile::LeftCutoff(tracepos t)
  297.     {
  298.     leftCutoff = t;
  299.     for(uint i=A;i<NUM_BASES;i++) trace[i]->LeftCutoff(t);
  300.     CalculateStats();
  301.     }
  302.  
  303. inline
  304. tracepos
  305. CTraceFile::RightCutoff(void)
  306.     {
  307.     return rightCutoff;
  308.     }
  309.  
  310. inline
  311. void
  312. CTraceFile::RightCutoff(tracepos t)
  313.     {
  314.     rightCutoff = t;
  315.     for(uint i=A;i<NUM_BASES;i++) trace[i]->RightCutoff(t);
  316.     CalculateStats();
  317.     }
  318.  
  319. inline
  320. tracepos
  321. CTraceFile::PrimerPosition(void)
  322.     {
  323.     return primerPosition;
  324.     }
  325.  
  326. inline
  327. CPeakList&
  328. CTraceFile::Peaks(void)
  329.     {
  330.     return peaks;
  331.     }
  332.  
  333. inline
  334. size_t
  335. CTraceFile::NumPeaks(void)
  336.     {
  337.     return peaks.Size();
  338.     }
  339.  
  340. inline
  341. CTrace<trace_t>*
  342. CTraceFile::operator[](enum_t whichTrace)
  343.     {
  344.     return trace[whichTrace];
  345.     }
  346.  
  347. inline
  348. void
  349. CTraceFile::NumBases(size_t nb)
  350.     {
  351.     numBases = nb;
  352.     }
  353.  
  354. inline
  355. size_t
  356. CTraceFile::NumBases(void)
  357.     {
  358.     return numBases;
  359.     }
  360.  
  361. inline
  362. void
  363. CTraceFile::Sequence(char* seq)
  364.     {
  365.     sequence = seq;
  366.     }
  367.     
  368. inline
  369. char*
  370. CTraceFile::Sequence(void)
  371.     {
  372.     return sequence;
  373.     }
  374.  
  375. inline
  376. void
  377. CTraceFile::Comments(char* newComment)
  378.     {
  379.     comments = newComment;
  380.     }
  381.  
  382. inline
  383. char*
  384. CTraceFile::Comments(void)
  385.     {
  386.     return comments;
  387.     }
  388.  
  389. inline
  390. tracepos*
  391. CTraceFile::BasePositions(void)
  392.     {
  393.     return basePositions;
  394.     }
  395.  
  396. #endif
  397.